/*!
 * \file DataStruct.h
 * \date 2021/01/24 10:57
 *
 * \author wanyzh
 * Contact: yizhao_wan@126.com
 *
 * \brief Source file of the CMesh class for the geometry and finite element mesh.
 *
 * TODO: long description
 *
 * \note
*/
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>

#include "CMesh.h"

using namespace std;

 vector<CPoint> &CMesh::getPoints()
{
	return m_P;
}

 vector<CElement>&CMesh::getElements()
{
	return m_T;
}

 vector<CPoint> &CMesh::getPbTest()
{
	return m_Points;
}

 vector<CElement> &CMesh::getTbTest()
{
	return m_Elements;
}
 vector<CPoint> &CMesh::getPbTrial()
{
	return m_PbTrial;
}

 vector<CElement> &CMesh::getTbTrial()
{
	return m_TbTrial;
}

 vector<vector<boundary_edge>> &CMesh::getAllBoundEdges()
{
	return m_allBoundaryEdges;
}
 vector<vector<boundary_node>> &CMesh::getAllBoundNodes()
{
	return m_allBoundaryNodes;
}

bool CMesh::readMeshFromFile(string filename)
{
	string suffixStr = filename.substr(filename.find_last_of('.') + 1);//get the suffix of the mesh file
	if (suffixStr == "neu")
	{
		int node_number;
		int ele_number;

		ifstream f(filename, ios::in);;
		if (!f.is_open())
		{
			cout << "\n error: open mesh file error\n";
			return false;
		}
		else
		{
			try
			{
				int ele_type;
				string s = "";
				for (int i = 0; i < 6; ++i)//skip 6 line for the file header
				{
					getline(f, s);
				}
				int NGRPS, NBSETS, NDFCD, NDFVL;
				f >> node_number >> ele_number >> NGRPS >> NBSETS >> NDFCD >> NDFVL;

				for (int i = 0; i < 3; ++i)//As >> operator does not handle the '\r' in the last line
				{
					getline(f, s);
				}

				m_Points.resize(node_number);
				m_Elements.resize(ele_number);

				//	double t1;
				for (int i = 0; i < node_number; ++i)
				{
					int index;
					double x, y;
					f >> index >> x >> y;//read the node coordinates
					m_Points[i] = CPoint(x, y, 0);
				}
				for (int i = 0; i < 3; ++i)
				{
					getline(f, s);
				}

				for (int i = 0; i < ele_number; ++i)
				{
					int id; // no check is done to see that they are monotonously increasing
					int ntype, ndp;
					f >> id >> ntype >> ndp;
					ele_type = ELE_TYPE(ndp);
					m_Elements[i].set_ele_type(ELE_TYPE(ndp));//set the element type, it is also the basis type
					for (int j = 0; j < ndp; ++j)
					{
						int list;
						f >> list;
						m_Elements[i].append(list - 1);
					}
				}

				getline(f, s);
				getline(f, s);


				char str1[30];
				int nTempBoundNum;
				char charTemp[10];
				char charTepm1[10];
				int GroupNum = 0;
				int temp1, temp2, temp3;
				for (int i = 0; i < NGRPS; ++i)
				{
					f >> charTemp >> str1 >> charTepm1;
					string str = str1;
					if (str == "GROUP")
					{
						f >> str1 >> temp1 >> charTepm1 >> nTempBoundNum >> charTemp >> temp2 >> charTemp >> temp3;
						for (int i = 0; i < 3; ++i)
						{
							getline(f, s);
						}
						int p[10];
						for (int i = 0; i < nTempBoundNum / 10; ++i)
						{
							f >> p[0] >> p[1] >> p[2] >> p[3] >> p[4] >> p[5] >> p[6] >> p[7] >> p[8] >> p[9];

						}
						int nTempGroupRow = nTempBoundNum - nTempBoundNum % 10;
						for (int i = 0; i < nTempBoundNum % 10; ++i)
						{
							int temp;
							f >> temp;
						}
						for (int i = 0; i < 3; ++i)
						{
							getline(f, s);
						}
						GroupNum++;
					}
				}

				int ntemp = 0;
				int data_type;//0=node,1=element
				for (int i = 0; i < NBSETS; ++i)
				{
					f >> str1 >> data_type >> nTempBoundNum >> temp2 >> temp3;


					if (data_type == 0)//the boundary is recorded by nodes
					{
						vector<boundary_node> boundarynodes;
						for (int i = 0; i < nTempBoundNum; ++i)
						{
							boundary_node bound_node;
							int index;
							f >> index;
							bound_node.glo_node_index = index - 1;
							boundarynodes.push_back(bound_node);
						}
						m_allBoundaryNodes.push_back(boundarynodes);
					}
					else if (data_type == 1)//the boundary is recorded by element
					{
						vector<boundary_edge> boundaryedges;
						for (int i = 0; i < nTempBoundNum; ++i)
						{
							boundary_edge boun_edge;
							int index;
							f >> index >> temp1 >> temp2;//read element node on the boundary
							boun_edge.ele_index = index - 1;
							switch (temp2)
							{
							case 1:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[1]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[2]);
								}
								else if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[1]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[2]);
								}

								break;
							case 2:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[1]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[2]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[4]);
								}
								else if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[1]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[2]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[4]);
								}
								break;
							case 3:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[0]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[4]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[0]);
								}
								else if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[3]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[4]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[6]);
								}
								break;
							case 4://case 4 only exists when the element has four edges

								if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[3]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[0]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (m_Elements[boun_edge.ele_index].getIds()[6]);
									boun_edge.p.second = (m_Elements[boun_edge.ele_index].getIds()[0]);
								}
								break;
							default:
								break;
							}
							boundaryedges.push_back(boun_edge);
						}
						m_allBoundaryEdges.push_back(boundaryedges);
					}
					ntemp += 1;
					for (int i = 0; i < 3; ++i)
					{
						getline(f, s);
					}
				}

				f.close();
				

				return true;
			}
			catch (exception& e)
			{
				cout << e.what() << endl;
				return false;
			}


		}
	}
	else if (suffixStr == "msh")
	{

	}
	else
	{
		cout << "\nUnknown mesh file" << endl;
		return false;
	}
	return false;

}



void CMesh::assignBoundNodesValue(vector<boundary_node>& bound, int index_vec,BOUNDARY_TYPE type, double(*value)(CPoint,double t))
{
	for (unsigned int i = 0; i < bound.size(); ++i)
	{
		bound[i].type = type;
		bound[i].value[index_vec]=(value);
		bound[i].mark.push_back(index_vec);
	}
}

void CMesh::assignBoundEdgesValue(vector<boundary_edge>& bound, int index_vec,BOUNDARY_TYPE type, double(*value)(CPoint,double t),double(*r)(CPoint,double t))
{
	if (Neumann==type)
	{
		for (unsigned int i = 0; i < bound.size(); ++i)
		{
			bound[i].type = Neumann;
			bound[i].value[index_vec]=value;//the value
			bound[i].mark.push_back(index_vec);
		}
	}
	else if (Robin==type)	
	{
		for (unsigned int i = 0; i < bound.size(); ++i)
		{
			bound[i].type = Robin;
			bound[i].value[index_vec]=(value);//the value
			bound[i].r[index_vec]=(r);//the r value
			bound[i].mark.push_back(index_vec);
		}
	}
	
}


void CMesh::generatePT()
{

	//generate P and T from Pb_test and Tb_test
	//Pb_test and Tb_test are read from gambit .neu file


	int index_p = -1;//global index for the point
	for (unsigned int i=0;i<m_Elements.size();++i)
	{
		if (ELE_3== m_Elements[i].get_ele_type())
		{			
			CElement T;
			for (int j = 0; j < m_Elements[i].get_ele_node_num(); j = j+1)
			{
				//find if the point m_Points[m_Elements[i].getIds()[j]] exist in P vector and return the index of the point 
				int index = getVectorIndex(m_P, m_Points[m_Elements[i].getIds()[j]]);

				//if the index==-1, then the point is not in P vector
				//and we add the point to P vector, and the global index of P is added one
				//the added ID for T is the new global index
				if (-1==index)
				{
					m_P.push_back(m_Points[m_Elements[i].getIds()[j]]);
					index_p += 1;
					T.append(index_p);
				}
				//if the point exist in the P vector
				//do not add the point in P vector
				//and the added ID for T is the index of the existed point in P vector
				else
					T.append(index);
				
			}
			T.set_ele_type(ELE_3);
			m_T.push_back(T);
		}

		else if ( ELE_4 == m_Elements[i].get_ele_type())
		{
			
			CElement T;
			for (int j = 0; j < m_Elements[i].get_ele_node_num(); j = j + 1)
			{
				int index = getVectorIndex(m_P, m_Points[m_Elements[i].getIds()[j]]);
				if (-1 == index)
				{
					m_P.push_back(m_Points[m_Elements[i].getIds()[j]]);
					index_p += 1;
					T.append(index_p);
				}
				else
					T.append(index);
			}
			T.set_ele_type(ELE_4);
			m_T.push_back(T);
		}
		
		else if (ELE_6 == m_Elements[i].get_ele_type())
		{
			
			CElement T;
			for (int j = 0; j < m_Elements[i].get_ele_node_num(); j = j + 2)
			{
				int index = getVectorIndex(m_P, m_Points[m_Elements[i].getIds()[j]]);
				if (-1 == index)
				{
					m_P.push_back(m_Points[m_Elements[i].getIds()[j]]);
					index_p += 1;
					T.append(index_p);
				}
				else
					T.append(index);
			}
			T.set_ele_type(ELE_3);
			m_T.push_back(T);
		}
		else if (ELE_9 == m_Elements[i].get_ele_type())
		{
			
			CElement T;
			for (int j = 0; j < m_Elements[i].get_ele_node_num() - 1; j = j + 2)
			{
				int index = getVectorIndex(m_P, m_Points[m_Elements[i].getIds()[j]]);
				if (-1 == index)
				{
					m_P.push_back(m_Points[m_Elements[i].getIds()[j]]);
					index_p += 1;
					T.append(index_p);
				}
				else
					T.append(index);
			}
			T.set_ele_type(ELE_4);
			m_T.push_back(T);
		}
	}
	//Pb_trial and Tb_trial
	m_PbTrial = m_Points;
	m_TbTrial = m_Elements;

	//exportPT2Vtk("mesh.vtk", m_P, m_T);
	
}
template<class T>
int CMesh::getVectorIndex(vector<T>  &v,T K)
{
	auto it = find(v.begin(), v.end(), K);

	// If element was found
	if (it != v.end())
	{

		// calculating the index
		// of K
		int index = it - v.begin();
		return index;
	}
	else {
		// If the element is not
		// present in the vector
		return -1;
	}
}
void CMesh::exportPT2Vtk(string filename,vector<CPoint> &P, vector<CElement> &T)
{
	map<ELE_TYPE, int> type_ele_vtk = {
		{ELE_3,5},
		{ELE_4,9},
		{ELE_6,22},
		{ELE_9,28}
	};
	ofstream f(filename);
	f << "# vtk DataFile Version 3.0" << endl;
	f << "Result" << endl;
	f << "ASCII" << endl;
	f << "DATASET UNSTRUCTURED_GRID" << endl;

	f << "POINTS \t" << P.size() << "\tdouble" << endl;
	for (unsigned int i = 0; i < P.size(); ++i)
	{
		f << P[i].getX() << "\t" << P[i].getY() << "\t" << P[i].getZ() << endl;
	}
	f << "CELLS \t " << T.size() << "\t" << (T[0].getIds().size() + 1)*  T.size() << endl;
	for (unsigned int i = 0; i < T.size(); ++i)
	{

		f << T[i].getIds().size() << "\t";

		if (T[i].get_ele_type() == ELE_3 || T[i].get_ele_type() == ELE_4)
		{
			for (unsigned int j = 0; j < T[i].getIds().size(); ++j)
			{
				f << T[i].getIds()[j] << "\t";
			}
		}
		else if (T[i].get_ele_type() == ELE_6)
		{
			for (unsigned int j = 0; j < T[i].getIds().size(); j += 2)
			{
				f << T[i].getIds()[j] << "\t";
			}
			for (unsigned int j = 1; j < T[i].getIds().size(); j += 2)
			{
				f << T[i].getIds()[j] << "\t";
			}
		}
		else if (T[i].get_ele_type() == ELE_9)
		{
			for (unsigned int j = 0; j < T[i].getIds().size() - 1; j += 2)
			{
				f << T[i].getIds()[j] << "\t";
			}
			for (unsigned int j = 1; j < T[i].getIds().size() - 1; j += 2)
			{
				f << T[i].getIds()[j] << "\t";
			}
			f << T[i].getIds()[8] << "\t";
		}
		f << endl;

	}
	f << "CELL_TYPES\t" << T.size() << endl;

	for (unsigned int i = 0; i < T.size(); ++i)
	{
		f << type_ele_vtk[T[i].get_ele_type()] << endl;
	}
}
